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We theoretically describe the charge ordering (CO) metal-insulator transition based on a quasi-one- 
dimensional extended Hubbard model, and investigate the finite temperature (T) properties across the 
transition temperature. Too- In order to calculate T dependence of physical quantities such as the spin 
susceptibility and the electrical resistivity, both above and below Too, a theoretical scheme is developed 
which combines analytical methods with numerical calculations. We take advantage of the renormalization 
group equations derived from the effective bosonized Hamiltonian, where Lanczos exact diagonalization 
data are chosen as initial parameters, while the CO order parameter at finite-T is determined by quantum 
Monte Carlo simulations. The results show that the spin susceptibility does not show a steep singularity 
at Tco, and it slightly increases compared to the case without CO because of the suppression of the spin 
velocity. In contrast, the resistivity exhibits a sudden increase at Too, below which a characteristic T 
dependence is observed. We also compare our results with experiments on molecular conductors as well 
as transition metal oxides showing CO. 

KEYWORDS: charge ordenng, molecular conductors, extended Hubbard model, exact diagonalization, 
bosonization, renormalization group, quantum Monte Carlo, transition metal oxides 



1. Introduction 

Charge ordering (CO) phase transition is now found ubiq- 
uitously in strongly correlated electron systems such as tran- 
sition metal oxides'' and molecular conductors.^' Its intuitive 
picture is simple: Electrons arrange spontaneously which re- 
sults in lowering of the symmetry from the underlying lat- 
tice structure, in order to gain repulsive Coulomb energy as 
in Wigner crystals. Nevertheless, the richness of this phe- 
nomenon is now widely recognized, owing to extensive ex- 
perimental as well as theoretical investigations in many types 
of compounds with different lattice geometries. 

From the theoretical point of view, in spite of numerous 
studies,-'' there remains a fundamental question not fully clar- 
ified yet, i.e., how the physical properties across the CO phase 
transition temperature (T), Tco. can be described. Typical 
mean-field analysis fails in reproducing finite-T phase tran- 
' sitions from a metallic state to paramagnetic insulating CO 
states, which are, however, often observed in many strongly 
correlated electronic materials. Therefore, treatments beyond 
the simple mean-field approximation, which consider the ef- 
fects of quantum fluctuation more properly, are necessary to 
describe such properties. 

In general, one-dimensional (ID) models can be treated by 
taking the quantum fluctuations into account in a controlled 
way, compared to higher dimensional systems, by numeri- 
cal as well as analytical methods. In fact, to describe CO 
in quarter-filled systems, the ID extended Hubbard model 
(EHM) including the repulsive Coulomb interactions of on- 
site, U, and intersite, V , has intensively been studied. Espe- 
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cially, its ground-state properties (see Fig. 1) are known in 
detail; a quantum phase transition occurs between the metal- 
lic Tomonaga-Luttinger (TL) liquid state characterized by the 
TL liquid parameter Kp, and the CO insulating state {Tqo = 
0).^*"^' However, the ID model does not show any phase tran- 
sition at finite T due to the enhanced low-dimensional fluctu- 
ations. 

On the other hand, the two-dimensional (2D) square lattice 
EHM at quarter-filling shows finite Tqo- Different techniques 
beyond the mean-field approximation have been applied to 
investigate the finite-T properties of this model, such as ex- 
act diagonalization (ED),^' slave-boson,^' dynamical mean- 
field,^'^' correlator projection,'"' and quantum Monte Carlo 
(QMC)"' methods. However, due to theoretical difficulties, 
the physical properties across Tqo such as the spin suscepti- 
bility and the electrical resistivity are not elucidated, and the 
interplay between spin and charge degrees of freedom is not 
fully explored yet. 

In this context, flie quasi-lD (QID) EHM, i.e., ID EHM 
chains coupled by the interchain Coulomb interaction V±, 
has recently been studied by analytical'^' '^' as well as nu- 
merical'"'' '^' methods by the present authors and co-workers, 
which shows a finite-T CO phase transition with concomitant 
metal-insulator transition at quarter-filling. In these studies, 
the interchain mean-field treatment'^' is applied and the re- 
sultant effective ID model is solved using different methods 
which properly take into account the large fluctuation effects: 
by the bosonization H- renormahzation group (RG) scheme in 
refs. 12 and 13, and by numerical techniques, i.e., the quan- 
tum transfer-matrix method in ref. 14 and the QMC method 
in ref. 15. 

In the former analytical approach, which has the advantage 
in investigating the critical regions, it is found that the V±- 
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Fig. 1. (Color online) Phase diagram of the quasi-one-dimensional ex- 
tended Hubbard model, obtained by the interchain mean- field theory. '2' In 
region (i), the CO state is stabilized when V±_ > with finite > 0. In 
region (ii), infinitesimal V±_ changes the TL Uquid state at Vs_ = into the 
CO state. In region (iii), where the CO state is obtained even in the purely 
ID case, the CO state is obtained at finite temperature by the interchain 
coupling. The boundary between (i) and (ii), and that between (ii) and (iii) 
are characterized by the value of Kp = 1/2 and that of Kp = 1/4 for the 
ID model. 

term considerably affects the stability of the CO state com- 
pared with that in the ID EHM.'^' Due to the dimensionality 
effect, Tco always becomes finite whenever the CO phase is 
realized, except for the critical point (line). Then, the ground 
state phase diagram of the ID EHM on the U-V plane is di- 
vided into three regions depending on the TL liquid parameter 
Kp, as in Fig. 1. They show different properties when is 
turned on: 

• Region (i) {Kp > 1/2] Finite value of is necessary 
to produce the CO state. 

• Region (ii) [1/4 < ifp < 1/2] Infinitesimal VI turns the 
system from TL liquid to CO insulator with finite Tco- 

• Region (iii) {Kp is not defined (CO insulating ground 
state for V±^ = 0)] Infinitesimal Vx makes Tco finite. 

Therefore, once V± is added, the CO phase critically enlarges 
from region (iii) to regions (ii)-i-(iii). In addition, the effects of 
the lattice dimerization along the chain direction and the frus- 
tration in the interchain interactions on the above QID model 
have been studied;'^' the lattice dimerization suppresses Tqq^ 
and the interchain frustration leads to a competition between 
CO states with different charge patterns.'^* 

In these studies, although the finite-T properties above Tqq 
are elucidated by taking advantage of the RG treatment which 
we will describe later, properties below Tco were hardly in- 
vestigated due to difficulties in determining the CO order pa- 
rameter in a self-consistent manner Meanwhile, the 'initial' 
values in the RG equations were taken as the bare TL param- 
eters derived from weak-coupling expansions of the EHM; 
such a procedure loses accuracy in general when the inter- 
action becomes large. Nevertheless, this drawback is not due 
to the phase Hamiltonian representation by the bosonization 
method nor the RG approach; it is because of the choice of 
the initial condition for the RG equations. 

On the other hand, the numerical quantum transfer matrix 
and QMC methods used in refs. 14 and 15 are suitable for the 
intermediate-to-strong coupling regime; by applying the in- 
terchain mean-field approach as well, the QID EHM and its 



extensions with different types of electron-lattice interaction 
are treated. Competitions and co-existences among different 
states, including the paramagnetic CO insulator, are clarified, 
and the T dependences of different order parameters as well 
as the charge and spin susceptibilities across the transition 
temperatures are computed. In particular, the QMC method 
could provide highly accurate results down to considerably 
low T. However, there are quantities which cannot be cal- 
culated by such numerical simulations directly, such as the 
electrical resistivity: one of the most basic information from 
experiments. 

Up to now, the above-mentioned analytical and numeri- 
cal approaches for the QID quarter-filled electron systems 
have been separately employed, being complimentary to each 
other In the present study, in contrast, we develop a com- 
bined method. The numerical ED data are implemented into 
the analytical bosonization -i- RG scheme, as initial conditions 
of the RG equations."^* As for the properties below Tco, the 
T dependence of the CO order parameter is calculated by the 
QMC method,'^* and then adopted to the scheme above. This 
combined method enables us to compute the T dependences 
of physical properties such as the spin susceptibility and the 
electrical resistivity across Tqo- 

The organization of this paper is as follows. In §2, our com- 
bined analytical and numerical method applied to the QID 
EHM is formulated. In §3, the T dependences of the physical 
quantities across the CO transition are shown. Section 4 is de- 
voted to the summary and discussions, including comparisons 
with experiments on different CO materials. Detailed descrip- 
tion of our theoretical approach is presented in Appendix A, 
and a benchmark check of this method applied to the ID Hub- 
bard model is given in Appendix B. 

2. Formulation 

In this section, the model and formulation for our calcula- 
tion are given. In § 2.1, first we explain the interchain mean- 
field approach to the QID EHM and the bosonization -i- RG 
method applied to the effective ID model, which were partly 
formulated in refs. 12 and 13. Then, in § 2.2 we explain how 
the numerical techniques are combined with this method. 

2.1 Bosonization + RG 

The QID EHM that we study consists of quarter-filled ID 
extended Hubbard chains coupled by the interchain interac- 
tion Vj_.'^^ The Hamiltonian iTgiD is given by 

i?QlD=^iT^D + ^-L, (1) 
3 

i,s 

+ Tiij^f Uij^l + V'^Uij rii+ij, (2) 

i i 

H± = V± X rii^j Ui^ji , (3) 

where and H± represent the ID EHM of the j-th chain 
and the interchain-coupling term, respectively. Here, t is the 
transfer integral between the nearest-neighbor sites along the 
chain direction; we do not consider the interchain hopping 
energy here. The operator c| j ^ creates an electron with spin 
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s =t or 4, at the i-th site in the j-th chain. The density op- 
erators are defined as riij^s = '^i,j,s^ij,s — 1/4 and riij = 
nij,t + nij,i. 

The interchain coupling term is treated within the inter- 
chain mean-field approximation'^^ as 

(4) 

Assuming the Wigner-crystal-type CO with two-fold period- 
icity, we postulate {riij) = {—ly^^n where n is the CO or- 
der parameter.'^' After this procedure, we obtain an effective 
ID Hamiltonian: 

i,s 

+ ^ rii^^ rii^l -\-V^^ni n^+i 

i i 

+ zV^nJ2i~'^yn^ + ^zLV^n^, (5) 

where the chain index is omitted. The number of the sites in 
the chain and that of the adjacent chains are expressed as L 
and z, respectively. 

By the bosonization method, the effective ID Hamiltonian 
eq. (5) can be expressed in terms of the bosonic fields. Then 
the Hamiltonian is separated into the charge part Hp and the 
spin part Ha as 



H=J dx{Hp + Ha) + \zLVi_n'', 
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where the phase variables satisfy [6'^(a;), (/)y(a;')] = 
msgu.{x — x')5ij,v with ji,u = pox a, and a is a short-distance 
cutoff. The parameters Kpo (iv'^o) and Vpo (Wao) are the bare 
TL-liquid parameter and velocity for the charge (spin) degree 
of freedom, respectively. They take non-universal values de- 
pending on the strength of interactions. 

In the charge sector, Hp, there appear two kinds of non- 
linear terms. The term proportional to cos Wp originates from 
the Sfcp umklapp scattering [fcF(= 7r/4) is the Fermi momen- 
tum], whose coupling constant 9^/4 is finite even in the purely 
ID EHM and it leads to the CO insulating ground state."*' ''-^"^ 
On the other hand, the cos 26 p term represents the 4fcF umk- 
lapp scattering process, whose coupling constant g^±^ is pro- 
portional to n when \g3i_\ <^ Trwpo-'^' This Akp umklapp 
scattering process is generated in the existence of CO, which 
can be understood by noting that the one-particle energy gap 
opens at ±2A:f because of the two-fold periodicity in the CO 
state and then the conduction band becomes effectively half 
filled. 

The spin part Ha is essentially the same as the effective 



Hamiltonian of the Heisenberg chain. The parameters Kao 
and g-i± in eq. (8) are not independent because of the spin- 
rotational SU(2) symmetry: 



KaO 



(9) 



This constraint still holds even under the scaling procedure. In 
the low-energy Umit, the gi±_ coupling is renormalized to zero 
and Kao reduces to unity. When the system has the SU(2) 
symmetry, it is known that physical quantities exhibit loga- 
rithmic (very slow) system-size or T dependences due to the 
presence of marginal operators.'^'' These characteristics can 
be captured by the RG analysis. 

The RG equations for the charge part Hp [eq. (7)] are given 

by 



-Kp{l) = -2GL(0 Klil) - 8Gl/,{l) KlU), 



(10) 



tGs^W = (2 - 2Kpil))G3^il) - G3±(Z)Gi/4(Z), (11) 



1 



-Gv4(0 = (2 - SKpil))G,/,il) - -GiAl)- 



(12) 



dl 
d_ 

We note that, in ref. 12, the above equations with Gs^ [1) = 
((?3_L = 0) were treated, corresponding to the absence of CO 
order parameter n, to investigate the instabiUty toward Tco- 
As for the spin part, the scaling equation for the coupling in 
Ha [eq. (8)] is 



and Ka{l) is determined as 



Ka{l) 



-2GL(0 



2GL(Z), 



1-Gi±(/)' 



(13) 



(14) 



following eq. (9). In the usual analysis, to obtain the pa- 
rameters {ifp(Z), G3_l(/), Gi/4(/)} for the charge part and 
G\_\_{1) for the spin part as the solutions of the RG equa- 
tions, the initial conditions are set from the bare parameter 
values as Kp{Qi) = Kpo, G3_l(0) = g3±/{TTVpo), Gi/4(0) = 
5(i/4/(27rfpo), and Gi_l(0) = gi±/{nVao)- Such initial val- 
ues can be calculated from the parameters of the original lat- 
tice model by considering the interaction processes near the 
Fermi level. In our previous studies on CO based on such 
procedure,^' '^-^^^ the third-order processes mediated by the 
states far from the Fermi level were crucial in deriving the 
Sfcp umklapp scattering gi/i-term, which triggers the CO in- 
sulating state. Note that the third-order virtual processes also 
play crucial roles for the spin degree of freedom."*' 

We also note that the RG equations above for the charge 
part are shown up to the one-loop level, while, on the other 
hand, that for the spin coupling Gi±{l) is shown up to the 
two-loop level, i.e., 0{G^j_). Since there are subtleties in 
deriving the two-loop RG equation based on the bosonized 
Hamiltonian, we follow the consideration given in ref. 23 and 
use the RG equation for Gij_ based on the Hamiltonian for 
the original fermion variables. 

When we apply the RG method to systems at finite T, 
the assumption of the scaUng invariance breaks down and 
the RG scaling is cut off at the scale I corresponding to 
the temperature T: I It = hi{Ct/T) with C be- 
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ing an 0(1) numerical constant. Then we can discuss the 
finite-T properties by terminating the scaling procedure at 
It, and use the values {Xp(It), G3±(^t), Gi/4(Zt)} and 
Gi_l(/t) as the T-dependent quantities; we write them as 
{Kp{T),Gzx{T),Gi,i{T)} and Gix{T) in the foUowing. 
These T-dependent parameters are set in the formulae for the 
physical quantities, as will be discussed in § 3. 

2.2 Numerical methods 

The procedure of setting the initial conditions of the RG 
equation mentioned above, i.e., to derive them based on the 
perturbation theory, generally looses accuracy in the strong- 
coupling region. For example, the ground state phase diagram 
of the ID EHM obtained in this way qualitatively agrees with 
the other numerical methods in the weak-coupling region, 
while the phase boundary between TL liquid and CO insulator 
deviates at strong-coupling.^' In the present study, instead, in 
order to obtain more accurate results even for stronger inter- 
actions, we make use of numerical data for finite size systems, 
as discussed in the following. Such an approach was recently 
proposed in ref . 1 8 to investigate the charge degree of freedom 
of the ID EHM, which reproduced the ground-state phase di- 
agram with high accuracy. Here we extend their method in 
order to calculate the finite-T properties by solving the RG 
equations in § 2.1 and stopping the scaling at It as mentioned 
above. 

For small systems with sites L, the TL parameters and 
can be computed by the Lanczos ED technique directly 
applied to the lattice Hamiltonian eq. (5). These can be used as 
initial values of the RG equations by taking advantage of the 
relation between the scaling variable I and the system size, I ~ 
In L (in the calculation we use Z = In L since the deviation is 
small for a different factor ). For the charge part, eqs. (10)- 

(12) , where the three parameters {Kp{l),Gi_\_{l),Giii{l)} 
are coupled in the form of differential equations, we can set 
Kp for available system sizes as initial values by fitting them 
to the RG equations to determine the solutions. As for the 
spin degree of freedom, the initial value for Gi±{l) in eq. 

(13) is determined from using the relation eq. (9). It will 
be discussed in § 3.1 that, for the calculation of the spin sus- 
ceptibility, we also make use of T-dependent spin velocity 
Va{T). However, the spin velocity cannot be determined re- 
Uably from the RG flow due to subtleties in deriving its RG 
equation. Therefore, instead, we use the standard polynomial 
finite size scaling for several system sizes L and extrapolate it 
to the size corresponding to In L ~ It, which provides Va (T). 
More detailed procedure for determining the T-dependent pa- 
rameters using the Lanczos ED data is given in Appendix A. 

We also combine numerical data for the CO order parame- 
ter n. It is finite for T < Tco, which affects the parameters 
of the system; especially ^sj, (Gsj,) becomes finite when n 
7^ 0. However, its T dependence cannot be obtained within 
the bosonization + RG scheme in a self-consistent manner, 
due to the ambiguity in the relationship between n and the 
phase variable.'^*' In order to overcome this difficulty, in the 
present study, n is numerically obtained by the QMC method 
in advance, and adopted to the above scheme. We employ 
the stochastic-series-expansion (SSE) method^^'^^' with the 
operator-loop update^^'^^' as a solver for the effective ID lat- 
tice Hamiltonian (5). At each T, the order parameter is itera- 



tively calculated by this QMC method until it converges.'^' 
We calculate systems with sizes up to T = 128 sites and 
checked that finite size effects are negligible down to low T 
that we discuss in this paper. The results are substituted into 
eq. (5) and treated as an 'external field', which reflect the ini- 
tial values of the RG equations. 

Summarizing, at each T, first the value of n is calculated 
by QMC, and then the initial conditions of the RG equations 
are collected using Lanczos ED by substituting the QMC data 
into the lattice Hamiltonian, and finally the RG equations (and 
the usual finite size scaling for the spin velocity) are solved by 
stopping the scaling at the scale I corresponding to the tem- 
perature T. This provides the finite-T values of the parameters 
to be input in the expressions of the physical quantities which 
will be described in the next section. 

3. Temperature Dependences of Spin Susceptibility and 
Electrical Resistivity 

In this section, the T dependences of the spin susceptibility 
Xcr(T) and the electrical resistivity p{T) are discussed. We 
focus on the regions (ii) and (iii) in Fig. 1 and study how the 
physical quantities behave for wide T ranges across Tco - For 
the intrachain parameter, we set {U /t, V/t) ~ (6.0, 2.5) for 
the region (ii) and (10.0, 4.0) for the region (iii). 

3.1 Uniform spin susceptibility 

The formula for the uniform spin susceptibihty Xct (T) in 
the RG scheme has been derived in refs. 29 and 30. The 
naive random-phase-approximation (RPA) gives Xct^^(^) — 
[2xo(T)/7rvF]/[l - Uxo{T)/{TrvF)\, where vp is the Fermi 
velocity and Xq^) is the spin susceptibility in the non- 
interacting case normaUzed as Xo(0) = 1. On the other hand, 
when the ID fluctuation effects are taken into account by the 
RG method, it is written as,^^-^"' 



XoiT) 



(15) 



TTVF 1 - [Gi±(T) + G4a{T)] xo(T) ■ 

Here, Gi_l(T) is the amphtude of the backward scattering, 
which is obtained by the RG scheme in § 2. The coupUng 
G4a-{T), on the other hand, represents the same-branch for- 
ward scattering in the spin channel.^^' It reflects the velocity 
of the spin excitations through the relation 



Gia{T) 



1 - 



Vf 



(16) 



then we can use the T-dependent spin velocity 7v(T) ob- 
tained by the finite size scaling introduced in § 2. We note 
that eqs. (15) and(16) reproduce the correct formula Xct(O) = 
2/[7rw<,(0)] in the T ^ limit. In refs. 29 and 30, the lin- 
earized dispersion was used in deriving the noninteracting 
susceptibility Xo (^) • Here, instead, in order to analyze x<t {T) 
in a wider T range, we use the dispersion of the tight-binding 
model to calculate Xo{T) (see also Appendix B). 

In Fig. 2, the results are shown for regions (ii) and (iii), 
where, in both cases, x<y{T) is enhanced below Tco, with- 
out any steep singularity at T = Tco- The main reason for 
the enhancement is as foUows. The CO order parameter n in- 
duces a gap formation at fc = ±2fcF = ±7r/2 in the energy 
dispersion, and hence the density of states at the Fermi energy 
is increased (the Fermi velocity is suppressed). This leads to 
the suppression of the spin velocity. 
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Fig. 2. (Color online) T dependence of the spin susceptibility Xa{T) 
(in the unit of l/{ta)), (a) for the region (ii) in Fig. \[{U/t,V/t) = 
(6.0, 2.5)] and (b) for the region (iii) [(U/t, V/t) = (10.0, 4.0)]. The 
solid curve shows the data for Vj_ = 0, without charge order at finite T. 



Fig. 3. (Color online) T dependence of the resistivity p(T) in arbitrary 
unit, (a) for the region (ii) and (b) for the region (iii) [same parameters as 
in Fig. 2]. The solid curve shows the data for V± = 0, without charge 
order at finite T. 



The same approach can be appHed to evaluate the spHtting 
of Knight shift, which is often observed in experiments to de- 
tect the CO transition. In the interchain mean-field approach, 
the Knight shift at the charge rich site and the poor site, 5+ 
and S-, are given by'^-* 



1 ± 



zV±n 



y/2t^ + {zV±ny 



(17) 



For small V±n/t, we obtain a simple relation S± oc 
Xcr{T)[l ± zV±n/2t]. Similarly, we can derive the formula 
of the nuclear relaxation rate, whose separation is also an ex- 
perimental evidence of CO: 



1 



iTi±)T 



1 ± 



zV±n 



^2t^ + {zV^nf 



(18) 



where the relaxation rate at the charge rich and poor 
sites are denoted by (Ti+)^^ and (Ti_)^^, respectively. 
Note that the discrepancy in the order of the factor 1 ± 
zVj_n/ \/2i? + (zV\nY' between (17) and (18) comes from 
the fact that the former expresses the local response under the 
uniform perturbation, whereas the latter is the local response 
under the local perturbation. Namely, Si cx Xcr(i, i'; 0) 
and (TiiT)^^ oc limc^^o IinXcr(«, where Si and T^,^ 

are the Knight shift and the NIVIR relaxation rate at the i-th 
site, and Xcr(*,i';'^) is the dynamical spin susceptibility in 
the site representation. 

3.2 Electrical resistivity 

Next we discuss the T dependence of the electrical resis- 
tivity p{T). Based on the memory function approach,^'"^^-* we 
perform the perturbation expansion with respect to g3± and 



in eq. (7) and then the formula for the resistivity is given 

by22) 



PiT) = 



1 



2a 



51/4 



2a 
6^ 



B^i^KpoAKpo) , (19) 



where B{x,y) — r(a;)r(y)/r(a; + y) is the beta function 
and = 'Vpo/i'^T) is the thermal coherence length. Here we 
note that the anomalous power-law behavior G'^^T'^^p~* and 
Gy^T^^^p^'^ can already be seen from the perturbationaly- 
obtained form, eq. (19). However this expression is valid 
only in the high-T region. In order to examine qualitative 
behaviors in the wide T range, we reinforce these expres- 
sions within the RG framework.^' By noting that the scal- 
ing dimension of the Akp and 8k p umklapp scatterings are 
2Kpo and 8X^0, the couplings are scaled (at tree level) 
as G3±(0 = G3±(0)cxp[(2-2Apo)Z] and G./^il) = 
Gi/4(0) exp [(2 — 8KpQ)l]. By inserting these relations into 
eq. (19), and replacing the bare parameter Kp^ by the renor- 
malized parameter Kp{l), we obtain the formula improved by 
the RG method. Then the T dependence can be computed as 



P{T) 



2tt^T 

VpO 



GlAT)B\Kp{T),Kp{T)) 



4Gl/^{T)B\4Kp{T)AKp{T)) 



(20) 



where the T-dependent coupUng constants are obtained by the 
procedure in § 2. 
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In Fig. 3, we show the results of p{T) thus calculated, for 
the regions (ii) and (iii). First, in the case of Vj_ =0 without 
CO at finite-T, the system shows a metallic behavior for the 
whole T range for region (ii), i.e., p{T) decreases with de- 
creasing T, whereas for region (iii), it is insulating below the 
T scale of the charge gap, i.e., p(T) increases with decreas- 
ing T. This behavior reflects the ground-state properties in the 
absence of the interchain coupUng, i.e., the TL liquid in the 
region (ii), while the CO insulating state in the region (iii). A 
noticeable point is that p{T) shows insulating behavior even 
without long range order of CO^^^ due to the low-dimensional 
fluctuation effect. 

When V±_ ^ 0, Tco becomes finite, then one can see a 
clear cusp in p{T) at T = Too- Just below Tco, p{T) shows 
a curve which is convex upward in the semi-log plot, reflect- 
ing the gap opening with the rapid growth of the CO order 
parameter n. At lower T, the curve turns convex downward, 
since at sufficiently low T, n becomes almost T independent 
and therefore the gap can be considered as a constant value, 
A; then an activation-type behavior p(T) cx exp (A/T) is 
expected. Such a behavior is common for all the parameters 
we have considered, as shown in Fig. 3. The abrupt change 
at T = Tco is indeed due to the emergence of 4A;F-umklapp 
scattering which originates from the gap in the energy disper- 
sion at ±2fcF owing to CO. This is in clear contrast with the 
behavior of the spin susceptibility Xct (T) shown in Fig. 2 with 
only a liny singularity at the transition. 

4. Summary and Discussions 

In the present paper, we have formulated a new theoretical 
framework to investigate the finite-T properties of QID elec- 
tron systems. The method has been applied to the CO phase 
transition in the QID EHM at quarter filling, where extended 
Hubbard chains are coupled via interchain Coulomb repulsion 
treated within the interchain mean-field approximation. 

In our scheme, we derive the bosonized Hamiltonian for 
the effective ID model and treat the RG equations, and by 
stopping the scaling procedure at the corresponding scale 
we obtain finite-T properties of the system. As for the ini- 
tial values of the RG equations we use the numerical results 
for the small systems obtained by the Lanczos ED method; 
these provide quantitatively good estimates even for strong 
coupling regime, in contrast with the conventional treatment 
where only the interaction processes between electrons near 
the Fermi energy are taken into account. In addition, QMC 
method is employed to calculate T dependence of the CO or- 
der parameter, which is necessary for the quantitative calcu- 
lations below Tco but difficult to determine in the bosoniza- 
tion + RG scheme. This framework enables us to calculate 
physical quantities across Tco such as the spin susceptibility 
Xu{T) as well as the electrical resistivity p{T) which is hard 
to calculate by the QMC simulation. 

The results show that, CO leads to an enhancement of 
X(t{T), mainly due to the reduction of the spin velocity, but 
without any steep singularity at T = Tco- Such features in 
the x<7 iT) curves shown in Fig. 2 are consistent with those 
calculated by purely numerical methods in refs. 14 and 15, 
which also showed a slight enhancement below Tco- In addi- 
tion, our results share similarities with the T = properties of 
the ID EHM, where x<t{T = 0) shows no singular behavior 
at the critical value of V for the emergence of CO, and contin- 



uously enhances when entering the CO phase. ^ On the other 
hand, as for p{T), the CO phase transition results in a sudden 
increase with a change of the slope with decreasing T, due 
to the generation of 4fcF-umklapp scattering which originates 
from the gap formation in the energy dispersion at ±2fcF . This 
is the first theoretical work, to the authors' knowledge, calcu- 
lating p{T) across the CO transition temperature starting from 
a microscopic correlated electronic model and taking full ac- 
count of the quantum and thermal fluctuation effects. 

Such results can be compared with the experiments, where 
XaiT) and p(T) are the most essential information for the 
bulk magnetic and electric properties of the system, measured 
most commonly. In fact, the sudden increase and the change 
of slope in p{T) are observed in a wide classes of compounds 
showing CO, even in materials which are apparently not ap- 
plicable to our QID model. On the other hand, there are dif- 
ferences in the behavior of Xa{T) at Tco among the com- 
pounds, as discussed below. Nevertheless, in many materials 
no noticeable change is seen in Xa{T), indicating a transition 
to the paramagnetic insulating CO phase, as in our calcula- 
tions. 

For example, in quarter-filled molecular conductors where 
many compounds show CO, QID materials such as 
(TMTTF)2X and (DCNQ1)2X (except the 7r-d mixed com- 
pound X =Cu) are candidates to be directly compared to 
our results, since ours can be considered as a microscopic 
model for their electronic properties. Several members of the 
(TMTTF)2X family showing CO, such as X=SbFG, AsFe, 
and Re04, indeed show kinks at Tco in their transport prop- 
erties.^^- A difference between our model and the situa- 
tion in the actual TMTTF compounds is that, the present 
calculation does not include the intrinsic lattice dimeriza- 
tion along the chains while it exists in the materials, which 
leads to another non-linear term in the bosonized Hamilto- 
nian.2^> In (DI-DCNQI)2Ag, even fliough recent studies^*- 
revealed that the system undergoes a more complex charge- 
lattice ordering than the simple CO we investigate in this pa- 
per, the T dependence in p{T) across T = Tco observed 
at high pressure resembles our calculated data, whereas the 
kink at T = Tco is smeared out at ambient pressure.-'^' All 
these QID molecular conductors show no anomaly in the bulk 
Xct(T) at T = Tco, while NMR measurements show the ap- 
pearance of atomic sites showing different Knight shifts and 
relaxation rates:"*^''*'^ these are also consistent with our anal- 
ysis. Recently, several QID compounds without dimerization 
have been synthesized where a CO transition is suggested, 
such as (o-DMTTF)2Br'*2) and (EDT-TTF-CONMe)2X \_X= 
AsFe and Br].''3.44) There, as in (DI-DCNQI)2Ag at ambient 
pressure, the anomaly at T = Tco in p(T) is not clearly seen 
possibly due to the strong fluctuation which may be underesti- 
mated in our calculation due to the interchain mean-field treat- 
ment. A noticeable point is that in (o-DMTTF)2Br,'*2^ Xa (T) 
shows a steep decrease at around T = Tco, distinct from 
the other QID materials above with a smooth variation there, 
whose origin remains unclear. 

Now there are many quasi-two-dimensional (Q2D) molec- 
ular conductors found to show CO. In such cases, the anomaly 
in p{T) and the characteristic curvature below T = Tco are 
again ubiquitously observed, while x<t{T) show different be- 
haviors from material to material. The latter diversity is due 
to the fact that, in the Q2D compounds, there exists a variety 
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in the anisotropy of transfer integrals originated from differ- 
ent molecular packings, which results in different anisotropic 
exchange couplings connecting the localized spins when CO 
is formed: For example, a spin gapped behavior is observed 
in a-(BEDT-TTF)2l3'^^' and ^"-(DODHT)2PF6'**'' due to the 
alternation in the exchange couplings along the charge rich 
sites. On the other hand, behavior analogous to our calcu- 
lation is seen in 6»-(BDT-TTP)2Cu(NCS)2,'^^''^^^ where the 
p(T) curve shows a rapid increase at around Tqq = 250 K, 
where Xa-^) shows no anomaly and below which the sys- 
tem is paramagnetic. In 6'-(BEDT-TTF)2RbZn(SCN)4, where 
the molecular packing in the 2D plane is the same, Xa {T) 
across Tco = 190 K is also paramagnetic; however, in 
p{T) a large jump is observed at T = Tco-"*'^ The dif- 
ference between the two ^-type compounds is in their CO 
pattern: the 'vertical stripe' in the former and the 'hori- 
zontal stripe' in the latter. This results in different man- 
ners in coupling to the lattice degree of freedom, leading 
to the second-order vs strong first-order nature of the CO 
phase transition. We note that most of the Q2D compounds 
show a second-order or a weak first-order phase transition, 
and 6'-(BEDT-TTF)2RbZn(SCN)4 is rather exceptional. An- 
other point we note here is that under pressure, in several 
Q2D compounds such as a-(BEDT-TTF)2l3^'" and ^-{meso- 
DMBEDT-TTF)2PF6,^''^^> p(T) curve at low T points toward 
a finite value extrapolated to T = 0; it does not diverge as in 
our calculation. Such a behavior is observed near the border 
between CO and uniform metallic phase without CO, and sug- 
gests the existence of a CO metallic phase. This is not reahzed 
in our calculation, where the CO phase is always insulating, 
therefore can be interpreted as the dimensionality effect in the 
transfer integrals. In fact, the CO metalhc phase is suggested 
in calculations on the 2D EHM.^^ 

Many transition metal oxides, even with Q2D or three- 
dimensional structures, show CO when the number of carri- 
ers become a fraction of the lattice sites. In such cases, again 
the sharp kink structure in p{T) is widely observed, e.g., 
in Nickelates, Manganites, Vanadates, and Iron based com- 
pounds.^^ However, the measurements for Xcr(T) shows a va- 
riety in their behavior, which is due to the same origin as in 
the Q2D molecular materials mentioned above: the variety in 
the anisotropy of the exchange couplings. One recent exam- 
ple where X<y{T) is continuous at the CO transition with a 
slight increase below Tco = 130 K is ^-Nai/3V205.52) This 
compound is a QID compound but has a compUcated crystal 
structure with a filling factor of 1/12; nevertheless the behav- 
ior, in the p{T) measurement as well for different pressures,^'*^ 
resembles our calculations. The results of calculation that we 
obtain can be applied to many material systems. 
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Appendix A: Combined Approacli of Numerical Anal- 
ysis, Bosonization, and Renormalization 
Group Method 

In this appendix we explain how to determine the T- 
dependent parameters needed to evaluate the formulae of eqs. 
(15) and (20), by combining the numerical results for the fi- 
nite size system and the RG equations. 

The TL parameters, Kp,K^, and the velocities, v^.v^, 
for the charge and spin degrees of freedom, respectively, in 
a finite L sites chain can be calculated exactly by the Lanc- 
zos ED technique using several standard relations. '^•''^^ The 
quantities above (we omit the superscript L) can be expressed 
as 

1/2 

Vo=\^\ , (A-la) 



1 



1/2 



tlJL 

TTK 



.1/2 



1/2 



(A- lb) 



where k and Xs are the compressibility and the spin suscepti- 
bility, respectively, whereas Dp and Dcr are the Drude weights 
for the charge and spin currents. For the finite size systems, k 
and Xs are given by 

K-i =^ [Eo{N^ + 1, AT^ + 1) + Eo{N^ + 1,N^- 1) 

-2EoiN^,N^)], (A-2a) 
=2L [Eo{N^ + 1, iV^ - 1) - £;o(iVt, N^)] , (A-2b) 

where Eq{N^, N^^) is the ground-state energy of the system 
with iV^ spin-up and spin-down electrons. For the estima- 
tion of the Drude weights, we use the relations 



TT d^Eoi^) 
L 



6=0 



d^En 



■'=0 

(A-3) 

where EQ{(j)) (Eo{(t>')) is the ground-state energy in the pres- 
ence of the charge (spin) flux $ = L(j) ($' = Lcp') through 
the ID ring. Then the TL liquid parameters and velocities are 
calculated exactly for finite size systems (up to 16 sites) using 
the Lanczos algorithm. 

To estimate the TL-liquid parameters at finite-T, we extend 
a recently-developed method combining the RG method with 
numerical results on finite size systems to study the charge 
degree of freedom in the ID EHM at quarter-filling. The 
essential effects of the nonlinear terms are to renormahze the 
parameters, namely, the TL-liquid parameters have expUcit 
size dependences. In ref. 18, it has been assumed that the L 
dependence of Kp is governed by the RG scaUng equations 
for the nonlinear term G1/4, and for several-system sizes 
L are used as the initial condition where the data is fitted to 
the scaling equations using the relation I = lnL. 

By following this idea, we determine the initial conditions 
of the RG equations {Kp{li),G3±{li),Gi/4^{li)} for the 
charge degree of freedom, eqs. (10)-(12), as follows, (i) We 
calculate k and Dp using eqs. (A-2) and (A-3) which give the 
TL-hquid parameter by eq. (A l), for three kinds of system 
size, K^=^, and K^='^^. (ii) We set Kp{h) = K^=^ 

and we tentatively substitute some values to G3_l(/i) and 
G'i/4(/i). By solving eqs. (10)-(12) with these tentative initial 
values, we calculate Kp(l2) and Kp{lz) where I2 = In 12, and 
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Fig. A l. (Color online) TL liquid parameter Ka of finite size systems 
L = 4,8,12,16 (open circles) for several choices of ((7/4, y/t).The 
solid lines denote the scaling trajectories of K^(X) with i = InL. 



^3 = In 16. (iii) We estimate {bKpf = {K p{h) ~ K^^^'^f ^ 
{Kp{ls) — Kp^^^)^. Through these steps, we optimize the 
initial values of G'3_l(^i) and 6*1/4(^1) so as to minimize the 
quantity SKp. Here we have assumed that G3±{li) (6*1/4(^1)) 
is odd (even) as a function of zV±n, by which the trivial limit 
that G3±{li) should reduce to zero for zV±n = is satisfied 
automatically. This assumption can be justified by consider- 
ing the perturbative forms. 

For the spin channel, the initial value of Gi± in eq. (13) 
is determined from the values of K„ of finite size systems. 
Owing to the spin-rotational SU(2) symmetry, the relation 
= [(l-hGi±)/(l-Gi_L)]^/2 (eq.(14)) still holds even for 
the finite size systems. By using this relation, Gi± for a finite 
size system is obtained from Ka, and can be used as the initial 
condition. This means that the initial condition of the nonlin- 
ear term for the spin channel can be determined without the 
fitting procedure and we can directly check the assumption in 
ref. 18. For the ID EHM, the TL liquid parameter in the 
finite size system {L — 4, 8, 12, 16) and the scaling trajectory 
given by eq. (13) with I =\nL are shown in Fig. A-1. The ini- 
tial value of Gi± is estimated with the least-square fit, i.e., in 
order to minimize {K^='^ - Ka{h)Y + (K^^^ - Ka{li)f + 
(^L=i2 _ Ka[l2)f + [K^^^^ - Ka{h)f. We can see that 
the scaling trajectory reproduces the numerical ED data very 
well, even for small system sizes. We note that the downward- 
convex dependence at large 1/L cannot be obtained by the 
one-loop RG, while the upward-convex 1/L dependence at 
small 1/L is nothing but the logarithmic singularity^'* and 
can be reproduced even in the one-loop level. 

The T dependence of the G^a (T) coupling can be obtained 
by using eq. (16), where we need the T dependence of the 
spin velocity. The RG equation for the spin velocity has not 
been obtained without ambiguity, since the correction is not 
logarithmic and the scaling invariance is not retained. Thus we 
utilize the simple size dependence of the velocity to examine 
its T dependence. The finite size scahng relation is expressed 
as 



Va{L) = Vcr{oo) + Cl/L + C2/L^ 



(A4) 



where ci and C2 are numerical constants. In order to obtain 
the T dependence of the spin velocity, we use eq. (A4) by 



substituting L 
constant. 



e ~ Ct/T with C being an (9(1) numerical 
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Fig. B-1. (Color online) T dependence of the uniform spin susceptibility 
Xa in unit of l/{ta) for the ID Hubbard model, for several choices of 
U/t. The filled circle and the open circle represent the results by "Lanc- 
zos+Bosonization+RG" and those by "Bosonization+RG", respectively 
(see text). These results are compared with QMC results (filled square for 
L = 64 and filled triangle for L = 128). 
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Fig. B-2. (Color online) T dependence of the normalized uniform suscep- 
tibility Xo(T) in the case of e{K) = —2t cos Ka (solid curve) and the 
susceptibility Xol{T) fov e{K) = vf{pK - kp) with \e{K)\ < Eo/2 
(dashed curve for Eq = 2t and dotted one for Eq = 4t). 



Appendix B: Spin Susceptibility in One-Dimensional 
Hubbard Model 

In this section, we demonstrate the validity and the ad- 
vantage of our combined analytical and numerical method 
by evaluating the uniform spin susceptibility of the ID Hub- 
bard model as a function of T. The present method (Lanc- 
zosH-BosonizationH-RG) is compared with the QMC method. 
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as well as the conventional method based on the RG the- 
ory (Bosonization+RG). We summarize the results in Fig. 
B-1, where the case for U/t = 2.0 (a), 3.0 (b), and 4.0 (c) 
are shown. Both the "Lanczos+Bosonization+RG" method 
and the usual "Bosonization+RG" method can reproduce the 
QMC results qualitatively. Especially, the present method 
gives more accurate results compared with the usual method 
even in the strong interaction cases. We note that, the dis- 
crepancy between the present results and the QMC method 
is mainly due to the formula of the spin susceptibility, eq. 
(15). The formula (15) is based on the RPA formalism with 
renormalized parameters, which overestimates the absolute 
value. ^''^ It should be noted that the noninteracting suscep- 
tibility Xo{T) is evaluated by using the cosine-band disper- 
sion even in the "Bosonization+RG" method instead of the 
linearized dispersion, in order to see the importance of choice 
of the initial values of the RG equations. In the case of the lin- 
earized dispersion, the noninteracting susceptibility is written 
as 

Eq 



Xol(T) = tanh 



AT' 



(Bl) 



where the energy dispersion is terminated as — £'o/2 < 
vpipK — kp) < Eq/2 with p = zt. This shows monotonic 
decrease with increasing T and the comparison between the 
cosine-band dispersion and the linearized dispersion is shown 
in Fig. B-2. 
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